package scu.maqiang.cfd;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES1L21;
import scu.maqiang.mesh.Mesh1L2;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class ConvectionDiffusion1D {

	public static void main(String[] args) {
		Mesh1L2 mesh1D = new Mesh1L2();
		mesh1D.line(0.0, 1.0, 10);
		double velX = 1.0;
		double k = 0.003;
		double t = 0.0, time = 0.001, dt = 0.00000001;
		FES1L21 fs = new FES1L21(mesh1D);
		SRMatrix K1 = new SRMatrix(fs.GetNdof());
		SRMatrix K2 = new SRMatrix(fs.GetNdof());
		SRMatrix C = new SRMatrix(fs.GetNdof());
		SRMatrix M = new SRMatrix(fs.GetNdof());
		fs.assembleConvection(new double[] {velX * dt}, BVPType.COMMON, C);
		double[] du = new double[fs.GetNdof()];
		fs.assembleMass(new double[] {1.0}, BVPType.LUMP_MASS, M);
		SRMatrix M0 = new SRMatrix(M);
		double[] u = new double[fs.GetNdof()];
		u[fs.GetNdof() - 1] = 1.0;
		double[] Mu = M0.mul(u);
		double[] Cu = C.mul(u);
		double[] K1u = K1.mul(u);
		double[] K2u = K2.mul(u);
		double[] RHS = new double[fs.GetNdof()];
		fs.applyBC_MBN(M, Direct.All, 1);
		fs.applyBC_MBN(M, Direct.All, 2);
		mesh1D.toTecplot("ConvectionDiffusion1D.dat", t, u);
		NewIterSSolver solver = new NewIterSSolver(M);
		while(t < time) {
			t += dt;
			System.out.println("time:" + t);
			Arrays.setAll(RHS, i -> Mu[i]-Cu[i]-K1u[i]-K2u[i]);
			fs.applyBC_MBN(RHS, Direct.All, 0.0, 1);
			fs.applyBC_MBN(RHS, Direct.All, 1.0, 2);
			solver.PCGSSOR(RHS, u, 1.5, 1);
			M0.mul(u, Mu);
			C.mul(u, Cu);
			K1.mul(u, K1u);
			K2.mul(u, K2u);
			mesh1D.toTecplot("ConvectionDiffusion1D.dat", t, u);
		}
	}

}
